For whatever reason I accidentally included the commit from my previous PR...

This is looking great!!

1. Regarding the dtypes, I think it should be done like numpy. I don't know how numpy does it haha but there must be somewhere in the documentation that says it.
2. Related to dtypes, maybe an .astype() or similar method should exist to help with these arithmetics?
3. Please also add a r1 / r2 operator if possible!
4. What happens if one raster has an incompatible amount of bands?

This is looking great!!

1. Regarding the dtypes, I think it should be done like numpy. I don't know how numpy does it haha but there must be somewhere in the documentation that says it.

Well, this is already following numpy's way. Under the hood it's simply an operation between numpy arrays so the output type of the array is set by numpy and then the Raster dtype is set from the numpy array. In my example with the numpy array for example, the output Raster has a type float64.

2. Related to dtypes, maybe an `.astype()` or similar method should exist to help with these arithmetics?

Mmmh, interesting idea! We sort of have this functionality already with the 'do-it-all' reproject function, which can convert dtype, except that right now if only the data type changes it won't do anything. So it could be a separate method yes.

3. Please also add a `r1 / r2` operator if possible!

I will! And r1 * r2.

4. What happens if one raster has an incompatible amount of bands?

If data.shape is not the same, an error is raised. So for now, both Rasters need to have the same amount of bands loaded into memory. We could also decide to use the first band by default, but it's so easy to load only a single or part of the bands with self.load that I would rather leave it as it is.

This is looking great!!

1. Regarding the dtypes, I think it should be done like numpy. I don't know how numpy does it haha but there must be somewhere in the documentation that says it.

Well, this is already following numpy's way. Under the hood it's simply an operation between numpy arrays so the output type of the array is set by numpy and then the Raster dtype is set from the numpy array. In my example with the numpy array for example, the output Raster has a type float64.

2. Related to dtypes, maybe an `.astype()` or similar method should exist to help with these arithmetics?

Mmmh, interesting idea! We sort of have this functionality already with the 'do-it-all' reproject function, which can convert dtype, except that right now if only the data type changes it won't do anything. So it could be a separate method yes.

3. Please also add a `r1 / r2` operator if possible!

I will! And r1 * r2.

4. What happens if one raster has an incompatible amount of bands?

If data.shape is not the same, an error is raised. So for now, both Rasters need to have the same amount of bands loaded into memory. We could also decide to use the first band by default, but it's so easy to load only a single or part of the bands with self.load that I would rather leave it as it is.

Sounds good with the dtype handling and band errors!

An .astype(whatever) method could just be added that runs .from_array with self._data.astype(whatever)

Actually, the question of dtypes is a problem for the subtraction...
Here's a minimal example that fails when subtracting a raster with dtypes float32 and uint8:

width = height = 5
transform = rio.transform.from_bounds(0, 0, 1, 1, width, height)
r1 = gu.Raster.from_array(np.random.randint(0, 255, (height, width)).astype('float32'), transform=transform, crs=None)
r2 = gu.Raster.from_array(np.random.randint(0, 255, (height, width), dtype='uint8'), transform=transform, crs=None)
r3 = r1 - r2
assert np.all(r3.data == r1.data - r2.data)

That's because running r1 - r2 will actually run r1 + -r2, and -r2 will be operated under the assumption of a uint8 type.

Why does it fail? Because -r2 doesn't work with a uint8 dtype? How does that work for uint8 ndarrays?

Why does it fail? Because -r2 doesn't work with a uint8 dtype? How does that work for uint8 ndarrays?

Because of this:
np.float32(10) - np.uint8(5) -> 5.0
np.float32(10) + -np.uint8(5) -> 261.0

It's actually a simple fix, I just need to make sure that the second argument has a compatible type with self, before doing the negation.

Nice!! Does * and / work or is this still in the works?

Nice!! Does * and / work or is this still in the works?

Still in the work. I'm working on the tests for the new methods right now.

@erikmannerfelt, if that version is ok, should we accept this PR and I'll create the multiplication/division in a new PR?